-
Notifications
You must be signed in to change notification settings - Fork 7
Add wind forcing and bottom drag tendency terms #219
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add wind forcing and bottom drag tendency terms #219
Conversation
|
Passes ctests on perlmutter and frontier, cpu and gpu. build and run sequence:
|
| | HTracersEdge | thickness-weighted tracers used for fluxes through edges. May be centered, upwinded or a combination of the two | ||
| | Del2TracersCell | laplacian of tracers on cells | ||
| | ZonalStressEdge | zonal component of wind stress on cells | ||
| | MeridStressEdge | meridional component of wind stress on cells |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think these should be ZonalStressCell and MeridStressCell, no? Later in the code these are defined on cells.
Also, I think having wind in the name would be good to differentiate from other potential sources of stress?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think that defining it on cells makes sense, given that eventually these fields are coupled on cell centers.
I think in one of our recent meetings we decided to keep the stress generic. It could be helpful to specify that it is at the top or surface though. Stress from both sea ice and wind currently aren't separable in MPAS-Ocean when they are passed from the coupler. We could decide to change this, of course.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@cbegeman I don't understand your second comment. It's not clear to me what external forced bottom stresses there may be (perhaps wave radiation stress is one)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I didn't mean to imply that we would have externally forced bottom stresses. Just that I think we're using these variables to refer to only stresses at the surface but the variable name quite generic.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Oh, I see. Thanks for the clarification. appending Surface does make sense for clarity. And I fully agree, these should be on cells, not edges.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think these should be ZonalStressCell and MeridStressCell, no? Later in the code these are defined on cells.
Yes, thanks for catching this.
To clarify how I came up with the names, in the code these variables are part of WindForcingAuxVars class. For me this means that the full name of ZonalStressCell is really WindForcingAuxVars.ZonalStressCell, and that makes it clear that this variable refers only to the wind stress (at the surface).
However, I see how this can be confusing, especially since the user docs do not mention the groupings of auxiliary variables. Moreover, there are no class prefixes in stream field names, so currently the name of ZonalStressCell used for I/O is WindStressZonal.
| MeshScalingDel4(Mesh->MeshScalingDel4), EdgeMask(Mesh->EdgeMask) {} | ||
|
|
||
| WindForcingOnEdge::WindForcingOnEdge(const HorzMesh *Mesh) | ||
| : SaltWaterDensity(1.026e3) {} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this meant to be the reference ocean density? If so, I think it should come from the namelist.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I changed this variable to come from the namelist
Tendencies:
Density0: 1026.0However, it seems to me that using the reference ocean density in this term is related to the Boussinesq approximation. We might want to change this.
| | Del2TracersCell | laplacian of tracers on cells | ||
| | ZonalStressEdge | zonal component of wind stress on cells | ||
| | MeridStressEdge | meridional component of wind stress on cells | ||
| | NormalStressEdge | normal component of wind stress on edge |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
generic comment - do we want this in AuxState? Would it make sense to have a 'forcing module' where this lives?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These are grouped in the WindForcingAuxVars class, so they aren't directly members of the aux state class. This is mentioned in the dev docs, but maybe it should also included be in the user docs.
| KOKKOS_FUNCTION void computeVarsOnEdge(int IEdge, int KChunk) const { | ||
| if (KChunk == 0) { | ||
| const int JCell0 = CellsOnEdge(IEdge, 0); | ||
| const int JCell1 = CellsOnEdge(IEdge, 1); | ||
| const Real ZonalStressEdge = | ||
| 0.5_Real * (ZonalStressCell(JCell0) + ZonalStressCell(JCell1)); | ||
| const Real MeridStressEdge = | ||
| 0.5_Real * (MeridStressCell(JCell0) + MeridStressCell(JCell1)); | ||
|
|
||
| NormalStressEdge(IEdge) = | ||
| Kokkos::cos(AngleEdge(IEdge)) * ZonalStressEdge + | ||
| Kokkos::sin(AngleEdge(IEdge)) * MeridStressEdge; | ||
| } | ||
| } |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should this be a more general edge reconstruction routine located elsewhere?
Also, this version is the anisotropic one we've moved away from for MPAS-Ocean. See E3SM-Project#6917.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I implemented the isotropic version and created a functor that can do both variants in HorzOperators.h. I added a namelist option to choose:
WindStress:
InterpType: Isotropic
| // Compute bottom drag | ||
| if (LocBottomDrag.Enabled) { | ||
| parallelFor( | ||
| {NEdgesAll, NChunks}, KOKKOS_LAMBDA(int IEdge, int KChunk) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
not sure why we need a loop over NChunks. I guess if Omega-0 is planning to do wind with multiple layers that's fine, but we'll have to change it later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I changed the bottom drag term to loop only over the edges in the bottom layer. For the wind stress term I kept the loop over multiple layers, since that's what is done in MPAS-O with some coefficient to distribute the stress that depends on the vertical coordinate. Once we have a vertical coordinate we can do the same in Omega.
|
Tested using barotropic gyre in polaris, which has an exact solution. Here I am using the same initial condition file for both. Both have the same viscosity, RK4 time step, same dt. MPAS-Ocean is running SW equations, with all vertical terms turned off. This is on a regular hex plane with 4200 cells. Omega Year 4 snapshot I would say that this is resounding support that this PR is working properly. We can retest if the stencil is altered per @cbegeman's comment above. The stencil difference from MPAS-O might cause the slight differences in the plots. I had to alter items by hand in I used gnu on perlmutter cpu, single node, for both. My run directories are here: Notes on how I ran this test:
|
|
@mark-petersen is the increase in error along the boundary in omega resulting from the 'partial slip' boundary in mpas-ocean? I can't recall if the boundary conditions are identical in mpas vs. omega |
@mark-petersen, I know you'll be away but would you have time to implement these changes in Polaris when you're back? |
I posted a comment along these lines on the relevant draft PR: E3SM-Project/polaris#286 (comment) |
|
Thanks @brian-oneill for adding the correction to the boundary masking, which is here: https://github.com/brian-oneill/E3SM/tree/omega/wind-forcing-bdry-edits. Now Omega and MPAS-O wind-forced results are much closer. Barotropic Gyre Global wind-forced case |
mark-petersen
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Based on the successful results above, I think this is now complete.
|
The plan is to merge this in, and then @brian-oneill will add the masking alterations in another PR to keep things clear. (see https://github.com/brian-oneill/E3SM/tree/omega/wind-forcing-bdry-edits) |
cbegeman
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Approving based on relatively small differences from MPAS-Ocean results and reasonable agreement with the analytic solution for a free slip boundary condition (expected increased drag associated with partial slip). Thanks, @mwarusz for your work on the code, @brian-oneill for prioritizing the masking, and @mark-petersen for testing.
katsmith133
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Approving based upon visual inspection and comparisons made by @mark-petersen between MPAS-O and Omega.
vanroekel
left a comment
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Approving based on visual inspection of code and testing from @mark-petersen, differences are small enough to proceed. Thanks @mwarusz and @brian-oneill!
|
@mark-petersen I fixed the errors and conflicts from the update to latest develop. Passes the unit tests on Chrysalis - all the changes were associated with the changes to Config interfaces and shouldn't impact anything else. |
|
Retested head. Passed all tests on Perlmutter CPU This error did not occur on the head of develop. My call sequence is as follows (see Frontier cpu section) build and run sequence:
|
|
@mark-petersen , I was able to reproduce the failure with the TIMESTEPPER test on Frontier using the |
|
Here's the issue: In |
5032675 to
45f937c
Compare
|
With previous addition from @brian-oneill: now passes all tests on Perlmutter CPU |












Adds wind forcing and bottom drag tendency terms. Includes some refactoring of
OceanTestCommon.hto make it easier to test 1D fields.This PR also partially addresses #178 because it removes the
Defaultsuffix from all auxiliary variables. However, I did not renameSshCellto justSsh. @xylarChecklist